C     Last change:  SWM  22 JUL 2002     4:35 pm
      SUBROUTINE LMG1AL(ISUM,ISUMI,LCA,LCIA,LCJA,LCU1,LCFRHS,LCIG,
     1                  ISIZ1,ISIZ2,ISIZ3,ISIZ4,ICG,NCOL,NROW,NLAY,
     2                  IN,IOUT,IFREFM)
C
C-----VERSION 1.0 25JUL2001 LMG1AL 
C
C     ******************************************************************
C     ALLOCATE STORAGE IN THE X ARRAY FOR LMG ARRAYS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
C
C1------PRINT A MESSAGE IDENTIFYING AMG PACKAGE
      WRITE(IOUT,500)
  500 FORMAT(1X,/1X,'LMG -- ALGEBRAIC MULTI-GRID SOLUTION PACKAGE',
     &     ', VERSION 1.1, 07/25/2001')
      
C
C2-------READ AND PRINT COMMENTS, STOR1, STOR2, STOR3, ICG
      CALL URDCOM(IN,IOUT,LINE)
      IF(IFREFM.EQ.0) THEN
        READ (LINE,505) STOR1, STOR2, STOR3, ICG
  505   FORMAT (3F10.0,I10)
      ELSE
        READ (LINE,*)  STOR1, STOR2, STOR3, ICG
      END IF
      WRITE(IOUT,510) STOR1, STOR2, STOR3
  510 FORMAT (1X,'STORAGE FACTOR FOR (A) AND (JA)= ',F6.3,/,
     &        1X,'STORAGE FACTOR FOR (IA), (U), AND (F)= ',F6.3,/,
     &        1X,'STORAGE FACTOR FOR (IG)= ',F6.3)
C
C2A-----IF CG CORRECTION NOT SET TO 0 OR 1, SET TO ZERO (NO CG)
      IF(ICG .NE. 0 .AND. ICG .NE. 1) ICG = 0
C
C3------ALLOCATE SPACE FOR THE AMG ARRAYS
C3------DOBULE PRECISION VARIABLES ARE: A, U, FRHS (LCA, LCU1, LCFRHS)
C3------INTEGER VARIABLES ARE: IA, JA, IG
      ISOLD=ISUM
      ISOLDI=ISUMI
      NRC=NROW*NCOL
      NRL=NROW*NLAY
      NCL=NCOL*NLAY
      NODES=NRC*NLAY
      NNA=NODES+2*(NCOL-1)*NRL+2*(NROW-1)*NCL+2*(NLAY-1)*NRC
      ISIZ1=STOR1*NNA+5*NODES
      ISIZ2=STOR2*NODES
      ISIZ3=STOR3*NODES
      ISIZ4=ISIZ2
      IF(ICG .EQ. 1)ISIZ4=ISIZ4+NODES
      LCA=ISUM                  
      ISUM=ISUM+ISIZ1  
      LCU1=ISUM                 
      ISUM=ISUM+ISIZ4  
      LCFRHS=ISUM               
      ISUM=ISUM+ISIZ4  
C3A-----INTERGERS
      LCIA=ISUMI
      ISUMI=ISUMI+ISIZ2
      LCJA=ISUMI
      ISUMI=ISUMI+ISIZ1
      LCIG=ISUMI
      ISUMI=ISUMI+ISIZ3
C
C4------CALCULATE AND PRINT THE SPACE USED IN THE Z AND IX ARRAYS
      IAMG=ISUM-ISOLD
      WRITE(IOUT,515) IAMG
  515 FORMAT(1X,I9,' ELEMENTS IN Z ARRAY ARE USED BY LMG')
      IAMG=ISUMI-ISOLDI
      WRITE(IOUT,520) IAMG
  520 FORMAT(1X,I9,' ELEMENTS IN IX ARRAY ARE USED BY LMG')
      WRITE(IOUT,525) ISIZ1, ISIZ4, ISIZ2, ISIZ3
  525 FORMAT(' SPACE ALLOCATED FOR [A] AND [JA]=',I9,/,
     &       ' SPACE ALLOCATED FOR [U] AND [FRHS]=',I9,/,
     &       ' SPACE ALLOCATED FOR [IA]=',I9,/,
     &       ' SPACE ALLOCATED FOR [IG]=',I9)
      RETURN
      END
C
C***********************************************************************
      SUBROUTINE LMG1RP(IN,MXITER,MXCYC,BCLOSE,DAMP,IOUTAMG,IOUT,
     &                  IFREFM,ICG,IADAMP,DUP,DLOW,HCLOSE) 
C
C-----VERSION 21FEB2001 LMG1RP
C     ******************************************************************
C     READ DATA FOR AMG
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
C     ------------------------------------------------------------------
C1-------INITIALIZE HCLOSE SO IT ISN'T UNDERFINED---erb 8/18/00
      HCLOSE = 0.0
      IADAMP = 0
      DUP=0.
      DLOW=0.
C
C2------READ MXITER,MXCYC,BCLOSE, DAMP, AND IOUTAMG
      IF(IFREFM.EQ.0) THEN
         READ (IN,500) MXITER,MXCYC,BCLOSE,DAMP,IOUTAMG
  500    FORMAT (2I10,2F10.0,I10)
      ELSE
         READ (IN,*)  MXITER,MXCYC,BCLOSE,DAMP,IOUTAMG
      END IF
C2------CHECK VALUE OF DAMP AND SEE IF ADAPTIVE DAMPING IS SELECTED.
C2------USE 1.E-6 AS A TOLERANCE WHEN CHECKING IF DAMP = -1 OR -2
C2------SET DAMP TO 1.0 FOR ALL OTHER VALUES, UNLESS THIS IS GREATER 
C2------THAN DUP, IN WHICH CASE SET DAMP TO DUP.
      IF(DAMP .LE. 0.)THEN
        IF(ABS(DAMP + 1.0) .LT. 1.E-6)THEN
          IADAMP = -1
        ELSEIF(ABS(DAMP + 2.0) .LT. 1.E-6)THEN
          DAMP=1.0
          IADAMP = -2
          READ(IN,*) DUP,DLOW
           IF(DUP .LT. 1.0)DAMP=DUP
        ELSE
          DAMP=1.0
        ENDIF
      ENDIF
C
C3------PRINT COMMENTS REGARDING SOLVER OPTIONS
      WRITE (IOUT,505)
  505   FORMAT (1X,///,36X,'SOLUTION BY THE ALGEBRAIC MULTIGRID METHOD',
     &        /,36X,42('-'))
      WRITE (IOUT,510) MXITER
  510 FORMAT (1X,20X,'MAXIMUM NUMBER OF CALLS TO AMG ROUTINE =',I9)
      WRITE (IOUT,515) MXCYC
  515 FORMAT (1X,25X,'MAXIMUM NUMBER OF CYCLES PER CALL =',I9)
      WRITE (IOUT,520) BCLOSE
  520 FORMAT (1X,10X,'SCALED L2 NORM OF RESIDUAL CRITERION FOR',
     &               ' CLOSURE =',E15.5)
      IF(ICG .EQ. 1)THEN
        WRITE(IOUT,521) ICG
      ELSE
        WRITE(IOUT,522) ICG
      ENDIF
  521 FORMAT (1X,30X,'CG-CORRECTION IS ACTIVE: ICG =',I9)        
  522 FORMAT (1X,26X,'CG-CORRECTION IS NOT ACTIVE: ICG =',I9)        
C3
      IF(IADAMP .EQ. -1)THEN
        WRITE(IOUT,525) IADAMP
      ELSEIF(IADAMP .EQ. -2)THEN
        WRITE(IOUT, 526) IADAMP
        WRITE(IOUT, 527) DUP
        WRITE(IOUT, 528) DLOW
      ELSE
        WRITE(IOUT,529) DAMP
      ENDIF
  525 FORMAT (1X,20X,'ADAPTIVE DAMPING USING COOLEY''S METHOD =', I9)
  526 FORMAT (1X,10X,'ADAPTIVE DAMPING USING RESIDUAL REDUCTION', 
     &               ' METHOD =',I9)
  527 FORMAT (1X,37X,'MAXIMUM VALUE OF DAMP =',E15.5)
  528 FORMAT (1X,37X,'MINIMUM VALUE OF DAMP =',E15.5)
  529 FORMAT (1X,41X,'DAMPING PARAMETER =',E15.5)
      IF(IOUTAMG .EQ. 0) WRITE(IOUT,530) IOUTAMG     
  530 FORMAT (1X,24X,'PRINTING FROM SOLVER IS SUPPRESSED =', I9)
      IF(IOUTAMG .EQ. 1) WRITE(IOUT,531) IOUTAMG     
  531 FORMAT (1X,28X,'PRINTING FROM SOLVER INCLUDES:',/,
     &          25X,'RESIDUALS BEFORE AND AFTER CYCLING =', I9)
      IF(IOUTAMG .EQ. 2) WRITE(IOUT,532) IOUTAMG     
  532 FORMAT (1X,28X,'PRINTING FROM SOLVER INCLUDES:',/,
     &           2X,'RESIDUALS BEFORE AND AFTER CYCLING, CP-TIMES,',
     &              ' AND STORAGE =',I9)
      IF(IOUTAMG .EQ. 3) WRITE(IOUT,533) IOUTAMG     
  533 FORMAT (1X,28X,'PRINTING FROM SOLVER INCLUDES:',/,
     &              'MESSAGES, RESIDUALS AFTER EACH CYCLE, CP-TIMES,',
     &              ' AND STORAGE =',I9)
C
C4------SET IOUTAMG TO A VALUE THAT AMG WANTS BY ADDING 10 TO THE VALUE.  
      IOUTAMG = 10+IOUTAMG
      RETURN
      END
C
C***********************************************************************
      SUBROUTINE LMG1AP(HNEW,IBOUND,CR,CC,CV,HCOF,RHS,A,IA,JA,U,FRHS,IG,
     &                  ISIZ1,ISIZ2,ISIZ3,ISIZ4,KITER,BCLOSE,DAMP,ICNVG,
     &                  KSTP,KPER,MXITER,MXCYC,NCOL,NROW,NLAY,NODES,
     &                  HNOFLO,IOUT,IOUTAMG,ICG,IADAMP,DUP,DLOW) 
C
C-----VERSION 1.1  25JUL2001 LMG1AP
C     ******************************************************************
C     THIS SUBROUTINE SETS UP THE MATRIX EQUATIONS IN A FORMAT THAT IS
C     COMPATIBLE WITH THE ALGEGRAIC MULTIGRID SOLVER PROVIDED BY GMD 
C     (AMG1R5).  THE 'AMG1R5' ALGORITHM IS THEN USED TO SOLVE FOR HEADS
C     OR SENSITIVITIES.
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*1 DIGIT(0:9)
      CHARACTER*2 PSUF
      CHARACTER*11 FNAME
      DOUBLEPRECISION HNEW,DZERO,RRHS,HHCOF,RSQ,FBAR,FMAX,DDAMP,DONE
      DOUBLEPRECISION Z,B,D,E,F,H,S
      DOUBLEPRECISION A(ISIZ1),U(ISIZ4),FRHS(ISIZ4),EPS,ECG1,ECG2,EWT2
      DOUBLEPRECISION ZHNEW, BHNEW, DHNEW, FHNEW, HHNEW, SHNEW
      PARAMETER (DZERO=0.D0, DONE=1.D0, IZERO=0)
      DIMENSION  IBOUND(NODES), HNEW(NODES),CR(NODES), CC(NODES),
     1   CV(NODES), HCOF(NODES), RHS(NODES), IA(ISIZ2),
     2   JA(ISIZ1),IG(ISIZ3)
      INCLUDE 'parallel.inc'
      DATA (DIGIT(I),I=0,9)/'0','1','2','3','4','5','6','7','8','9'/
C     ------------------------------------------------------------------
C
C1------INITIALIZE SCALAR VARIABLES VARIABLES AND CLEAR ARRAYS
C1------NCOUNT IS A COUNTER FOR THE NUMBER OF ACTIVE NODES
C1------NA AND NJ ARE INDEXES FOR THE MATRIX STORAGE USED BY AMG
      ICNVG=0
      RSQ=DZERO
      FMAX=DZERO
      FBAR=DZERO
      NCOUNT=0                  
      NRC=NROW*NCOL
      NA=1                      
      NJ=1                      
C1
      DO 100 I= 1, ISIZ1
        A(I)=DZERO
        JA(I)=IZERO
  100 CONTINUE
      IF(ICG .EQ. 0)THEN
        DO 105 I = 1, ISIZ2
          IA(I)=IZERO
          FRHS(I)=DZERO
          U(I)=DZERO
  105   CONTINUE
      ELSE
        DO 106 I = 1, ISIZ2
          IA(I)=IZERO
  106   CONTINUE
        DO 107 I = 1, ISIZ4
          FRHS(I)=DZERO
          U(I)=DZERO
  107   CONTINUE
      ENDIF
      DO 110 I=1,ISIZ3
        IG(I)=IZERO
  110 CONTINUE
C
C2------LOOP THROUGH ALL NODES IN THE GRID AND SET UP MATRIX EQUATIONS.
C2------THIS LOOP STRUCTURE AND INDEXING IS IDENTICAL TO THAT OF PCG2 
C2------AND IS BLATANTLY COPIED FROM HILL, 1990.
      DO 115 K=1,NLAY
      DO 115 I=1,NROW
      DO 115 J=1,NCOL
        N = J + (I-1)*NCOL + (K-1)*NRC
C2------CALCULATE 1 DIMENSIONAL SUBSCRIPT OF CURRENT CELL AND
C2------SKIP CALCULATIONS IF CELL IS INACTIVE
C2
C2------CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR LOCATING THE 6
C2------SURROUNDING CELLS
        NRN = N + NCOL
        NRL = N - NCOL
        NCN = N + 1
        NCL = N - 1
        NLN = N + NRC
        NLL = N - NRC
C2
C2------CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR CONDUCTANCE TO THE 6
C2------SURROUNDING CELLS.
        NCF = N
        NCD = N - 1
        NRB = N - NCOL
        NRH = N
        NLS = N
        NLZ = N - NRC
C2
C2------GET CONDUCTANCES TO NEIGHBORING CELLS
C2------NEIGHBOR IS 1 ROW BACK
        B = DZERO
        BHNEW = DZERO
        IF (I.NE.1) THEN
          B = CC(NRB)
          BHNEW = B*(HNEW(NRL)-HNEW(N))
        ENDIF
C2
C2------NEIGHBOR IS 1 ROW AHEAD
        H = DZERO
        HHNEW = DZERO
        IF (I.NE.NROW) THEN
          H = CC(NRH)
          HHNEW = H*(HNEW(NRN)-HNEW(N))
        ENDIF
C2
C2------NEIGHBOR IS 1 COLUMN BACK
        D = DZERO
        DHNEW = DZERO
        IF (J.NE.1) THEN
          D = CR(NCD)
          DHNEW = D*(HNEW(NCL)-HNEW(N))
        ENDIF
C2
C2------NEIGHBOR IS 1 COLUMN AHEAD
        F = DZERO
        FHNEW = DZERO
        IF (J.NE.NCOL) THEN
          F = CR(NCF)
          FHNEW = F*(HNEW(NCN)-HNEW(N))
        ENDIF
C2
C2------NEIGHBOR IS 1 LAYER BEHIND
        Z = DZERO
        ZHNEW = DZERO
        IF (K.NE.1) THEN
          Z = CV(NLZ)
          ZHNEW = Z*(HNEW(NLL)-HNEW(N))
        ENDIF
C2
C2------NEIGHBOR IS 1 LAYER AHEAD
        S = DZERO
        SHNEW = DZERO
        IF (K.NE.NLAY) THEN
          S = CV(NLS)
          SHNEW = S*(HNEW(NLN)-HNEW(N))
        ENDIF
C2    
C2------SKIP CALCULATIONS AND MAKE CELL INACTIVE IF ALL
C2------SURROUNDING CELLS ARE INACTIVE AND SET HNEW TO HNOFLO
        IF (B+H+D+F+Z+S.EQ.0. .AND. IBOUND(N) .NE. 0) THEN
          IBOUND(N) = 0
          HNEW(N) = HNOFLO
        ENDIF
C2
C2------CALCULATE L2 NORM FOR ACTIVE CELLS.  
        IF(IBOUND(N) .GT. 0)THEN
          RRHS = RHS(N)
          HHCOF = HNEW(N)*HCOF(N)
          RSQ = RSQ + (RRHS - ZHNEW - BHNEW - DHNEW - HHCOF - FHNEW - 
     &          HHNEW - SHNEW)**2
        ENDIF
C3
C3------SET THE DIAGONAL AND RHS AS A POSITIVE VALUE.  CHECK FOR ACTIVE 
C3------CELLS.  IF INACTIVE OR CONSTANT HEAD, SET DIAGONAL TO 1.0, OFF-
C3------DIAGONALS TO ZERO, AND ADJUST RHS ACCORDINGLY.  
C3------ACCUMULATE THE AVERAGE ABSOLUTE VALUE AND MAXIMUM VALUE OF THE 
C3------RHS VECTOR FOR ALL ACTIVE CELLS.  THIS IS USED TO SCALE THE 
C3------THE CLOSURE CRITERIA.
        E=Z+B+D+F+H+S-HCOF(N)
        IF(IBOUND(N) .GT. 0)THEN
          FRHS(N)=-RHS(N)
        ELSE
          E=DONE
          Z=DZERO
          B=DZERO
          D=DZERO
          F=DZERO
          H=DZERO
          S=DZERO
          FRHS(N)= HNEW(N)
        ENDIF
C3
        FMAX=MAX(FMAX,ABS(FRHS(N)))
        IF(IBOUND(N) .NE. 0)THEN
          FBAR=FBAR+ABS(FRHS(N))
          NCOUNT=NCOUNT+1
        ENDIF
C4
C4------USE A TYPE OF SKYLINE STORAGE WHICH AMG WANTS (SEE THE README 
C4------FILE DISTRIBUTED WITH THE AMG1R5 ALGORITHM). FOR EACH ROW, 
C4------START AT THE DIAGONAL.  THEN PICK UP ELEMENTS MOVING FROM LEFT 
C4------TO RIGHT (SKIPPING THE DIAGONAL). 'A' STORES ALL THE COEFFICENTS
C4------OF THE MATRIX IN A SINGLE VECTOR.  'IA' STORES THE LOCATION OF 
C4------THE DIAGONAL ELEMENTS WITHIN 'A', AND 'JA' STORES THE LOCATION 
C4------OF THE OFF-DIAGONAL ELEMENTS WITHIN 'A'
        A(NA) = E
        IA(N) = NA
        JA(NJ) = N
        NA = NA + 1
        NJ = NJ + 1
C4
        IF(K .NE. 1)THEN
          A(NA) = -Z
          JA(NJ) = N - NRC
          NA = NA + 1
          NJ = NJ + 1
        ENDIF
        IF(I .NE. 1)THEN
          A(NA) = -B
          JA(NJ) = N - NCOL
          NA = NA + 1
          NJ = NJ + 1
        ENDIF
        IF(J .NE. 1)THEN
          A(NA) = -D
          JA(NJ) = N - 1
          NA = NA + 1
          NJ = NJ + 1
        ENDIF
        IF(J .NE. NCOL)THEN
          A(NA) = -F
          JA(NJ)= N + 1
          NA = NA + 1
          NJ = NJ + 1
        ENDIF
        IF(I .NE. NROW)THEN
          A(NA) = -H        
          JA(NJ)= N + NCOL
          NA = NA + 1
          NJ = NJ + 1
        ENDIF
        IF(K .NE. NLAY)THEN
          A(NA) = -S
          JA(NJ)= N + NRC
          NA = NA + 1
          NJ = NJ + 1
        ENDIF
C4------STORE INITIAL GUESS
        U(N) = HNEW(N)     
  115 CONTINUE
C
C4
C4------SET THE FINAL LOCATIONS FOR THE COEFFICIENTS IN THE 'A' MATRIX
C4------REMEMBER THAT 'NA' AND 'NJ' HAVE ALWAYS BEEN INDEXED 1 AHEAD.
      IA(N+1) = NA              
      NA = NA - 1               
      NJ = NJ - 1               
C
C5
C5------SET THE DIMENSIONS USED FOR THESE VARIABLES AND PASS TO THE AMG
C5------SOLVER (AMG1R5).  IT USES THIS EXTRA SPACE TO WORK IN.
C5------ALSO SET VARIABLES FOR SOLVER CONTROL
      NDA = ISIZ1
      NDIA= ISIZ2
      NDJA= ISIZ1
      NNU= NODES
      NDU= ISIZ4
      NDF= ISIZ4
      NDIG= ISIZ3
C5
      MATRIX=22
      ISWTCH = 4
      IFIRST = 10
C
C6------CONTROL PRINTING FROM THE SOLVER.  IF OUTPUT IS LIMITED THEN
C6------OPEN TEMPORARY FILE ON FIRST AVAILABLE UNIT NUMBER. AMG MESSAGES
C6------ARE WRITTEN TO THIS FILE.  LATER, CLOSE AND DELETE THIS FILE.
C6------IF NO UNIT NUMBERS ARE AVAILABE, SEND AMG MESSAGES TO LIST FILE.  
      IPRINT = 10000+IOUT*100+IOUT
      IF(IOUTAMG .LE. 12)THEN
        FNAME='lmg_err.tmp'
        IF(NUMPROCS .GT. 1) THEN
          ITEN=MYID/10
          IONE=MOD(MYID,10)
          PSUF(1:1)=DIGIT(ITEN)
          PSUF(2:2)=DIGIT(IONE)
          FNAME='lmg_err.p'//PSUF
        ENDIF
        ILMGTMP=IGETUNIT(1,1000)
        IF(ILMGTMP .NE. -1)THEN
          OPEN(ILMGTMP, FILE=FNAME,STATUS='UNKNOWN')      
          IPRINT = 10000+IOUT*100+ILMGTMP
        ENDIF
      ENDIF 
C
C7------SET UP CLOSURE CRITERIA THEN CHECK CONVERGENCE 
      FBAR=FBAR/FLOAT(NCOUNT)
      EPS=BCLOSE/FMAX*FBAR
C
C7------IF CONVERGENCE CRITERIA IS MET, SET ICNVG=1 AND SKIP AMG SOLVER 
      IF(DSQRT(RSQ) .LE. EPS*FMAX .AND. KITER .GT. 1)THEN
        ICNVG = 1
        GOTO 150
      ENDIF
C7
C7------IF OUTPUT ISN'T SUPRRESSED, PRINT SCALING FACTORS FOR CLOSURE
      IF(IOUTAMG .NE. 10) WRITE(IOUT,500) FBAR,FMAX,EPS*FMAX
  500 FORMAT(//,3X,'AVERAGE ABSOLUTE VALUE OF THE RHS VECTOR=',E15.5, /,
     &       12X, 'MAXIMUM VALUE OF THE RHS VECTOR=',E15.5,/,
     &       10X, 'CLOSURE CRITERION FOR THE L2 NORM=',E15.5,/)
C
C8
C8------ADAPTED FROM AUXILLERY ROUTINE DISTRIBUTED WITH AMG1R5
C8------NOTE: MXCYC IS ADDED TO NCYC TO CONTROL MAX. NUMBER OF CYCLES
C8------AND 2000*ICG IS ADDED TO NCYC TO CONTROL CG-CORRECTION.
C8------FOR ZERO-PARAMETERS STANDARD VALUES ARE USED BY AMG1R5
      LEVELX = 0
      NCYC   = 10300 + 2000*ICG + MXCYC                 
      MADAPT = 0
      NRD    = 0
      NSOLCO = 2
      NRU    = 0
      ECG1   = 0.D0
      ECG2   = 0.25D0
      EWT2   = 0.35D0
      NWT    = 2
      NTR    = 0
C8
C8------CALL THE ACTUAL AMG SOLVER AS PROVIDED BY GMD.  RETURNS THE 
C8------SOLUTION IN VECTOR 'U'
      CALL AMG1R5(A,IA,JA,U,FRHS,IG,
     +            NDA,NDIA,NDJA,NDU,NDF,NDIG,NNU,MATRIX,
     +            ISWTCH,IOUTAMG,IPRINT,
     +            LEVELX,IFIRST,NCYC,EPS,MADAPT,NRD,NSOLCO,NRU,
     +            ECG1,ECG2,EWT2,NWT,NTR,
     +            IERR)
C
C9------CHECK FOR SOLVER ERRORS. IF THERE IS A SOLVER ERROR, REPORT IT
C9------AND STOP EXECUTION.
C9------IF THE SOLVE WAS SUCCESSFUL, THEN APPLY DAMPING AND RESTORE THE
C9------SOLUTION VECTOR IN HNEW.
C9------CLOSE AMG TEMPORARY FILE IF NECCESSARY
      IF(IERR .GT.0)THEN
        WRITE(IOUT,600) IERR
        IF(IERR .EQ. 1 .OR. IERR .EQ. 3) WRITE(IOUT,601)
        IF(IERR .EQ. 2 .OR. IERR .EQ. 4 .OR. IERR .EQ. 5)
     &       WRITE(IOUT,602)
        IF(IERR .EQ. 6) WRITE(IOUT,603)
        STOP
      ENDIF
  600 FORMAT(/,1X,'AMG SOLVER ERROR --> IERR= ',I3,/,
     &         1X,'SOLVE UNSUCCESSFUL.  CHECK VALUE OF IERR')
  601 FORMAT(/,1X,'STORAGE TO SMALL FOR [A] OR [JA]',/
     &         1X,'INCREASE SIZE OF [STOR1] IN LMG INPUT FILE')
  602 FORMAT(/,1X,'STORAGE TO SMALL FOR [IA], [U], OR [F]',/
     &         1X,'INCREASE SIZE OF [STOR2] IN LMG INPUT FILE')
  603 FORMAT(/,1X,'STORAGE TO SMALL FOR [IG]',/
     &         1X,'INCREASE SIZE OF [STOR3] IN LMG INPUT FILE')  
C9
C9------SET DAMPING.  IF USING AN ADAPTIVE DAMPING STRATEGY THEN
C9------IMPLEMENT DICK COOLEY'S DAMPING STRATEGY W/ HUYAKORN'S CHANGE
C9------OR IMPLEMENT RESIDUAL REDUCTION METHOD OF ADAPTIVE DAMPING
      DDAMP=DAMP
      IF(IADAMP .EQ. -1)THEN
        DHMAX=U(1)-HNEW(1)
        DO 120 N=2,NODES 
          TMP=U(N)-HNEW(N) 
          IF(ABS(DHMAX) .LT. ABS(TMP))DHMAX=TMP
  120   CONTINUE
        CALL ADAMP1(DDAMP,DHMAX,KITER)
C
      ELSEIF(IADAMP .EQ. -2)THEN
        CALL ADAMP2(RSQ,DDAMP,DAMP,DUP,DLOW,KITER)
      ENDIF
      DAMP=DDAMP
C9
      DO 130 N=1,NODES
        HNEW(N) = DDAMP*U(N) + (DONE-DDAMP)*HNEW(N)
  130 CONTINUE

C9
  150 IF(IOUTAMG .LE. 12 .AND. ILMGTMP .NE. -1) CLOSE(UNIT=ILMGTMP,
     &                                          STATUS='DELETE')    

      RETURN
      END
C
C***********************************************************************
      SUBROUTINE ADAMP1(WK,DHMAX,KITER)
C
C-----VERSION 1.0  31JAN2001 ADAMP1
C     ******************************************************************
C     THIS SUBROUTINE CALCULATES THE DAMPING PARAMETER USING DICK 
C     COOLEY'S SCHEME W/ HUYAKORN'S MODIFICATION
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DOUBLEPRECISION WK
      SAVE DHOLD

C1------INITIALIZE VARIABLES.  CALCULATE "S" IF NOT ON FIRST ITERATION
      S=1
      IF(KITER .GT. 1)S=DHMAX/WK/DHOLD
C2
C2------SET APPROPIATE DAMPING BASED ON "S".  "S" REPRESENTS THE CHANGE
C2------IN HEAD SCALED BY THE DAMPING.  NOTE THAT NEGATIVE VALUES 
C2------INDICATE THAT THE SIGN CHANGED FOR THE MAXIMUM HEAD CHANGE.
C2------FOR -1 < S < 0, THE MAX HEAD CHANGE CHANGED SIGN, BUT WAS 
C2------SMALLER IN MAGNITUDE FOR THE CURRENT ITERATION, SO DAMPING IS
C2------NOT THAT AGGRESSIVE (MINIMUM POSSIBLE IS 1/2).  
C2------FOR 0 < S < 1, THERE IS NO SIGN CHANGE AND THE MAXIMUM CHANGE  
C2------IS LESS FOR THIS ITERATION, WHICH IS AN INDICATOR OF GOOD 
C2------CONVERGENCE, SO NO DAMPING IS APPLIED (WK=1).  FOR S > 1, SAME
C2------IDEA, BUT THE MAX HEAD CHANGE WAS GREATER THIS ITERATION.
C2------FOR S < -1, MAX HEAD CHANGED SIGN, AND WAS GREATER IN MAGNITUDE
C2------WHICH INDICATES THE SOLUTION IS OSCILLATING AND DIVERGING, SO
C2------DAMPING IS AGGRESSIVE (MAXIMUM POSSIBLE IS 1/2).
      IF(S .GE. -1.)THEN
        WK=(3.+S)/(3.+ABS(S))
      ELSE
        WK=1./2./ABS(S)
      ENDIF
C3
C3-----STORE MAXIMUM HEAD CHANGE FOR NEXT ITERATION
      DHOLD=DHMAX
C     
      RETURN
      END
C
C***********************************************************************
      SUBROUTINE ADAMP2(RSQ1,DDAMP,DAMP,DUP,DLOW,KITER)
C
C-----VERSION 1.1  22JUL2002 ADAMP2
C     ******************************************************************
C     THIS SUBROUTINE CALCULATES THE DAMPING PARAMETER USING THE 
C     RESIDUAL REDUCTION METHOD
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      LOGICAL OSCIL
      DOUBLEPRECISION RSQ1,RSQ2,RSQ3,DDAMP
      DATA ICONST, IMULT, IPER /1283,106,6075/
      SAVE RSQ2,RSQ3,NRAN,DAMP3

C1
C1------CALCULATE RESIDUAL AND INITIALIZE VARIABLES FOR FIRST ITERATION
      RSQ1=DSQRT(RSQ1)
      OSCIL = .FALSE.
      IF(KITER .EQ. 1)THEN
        RSQ3=2.*RSQ1+1.
        DAMP3=DAMP
        NRAN=1
C2
C2------CALCULATE THE RELATIVE REDUCTION IN THE RESIDUAL FOR THIS
C2------ITERATION ("RRED") SCALED BY THE DAMPING.  
C2A-----IF THE REDUCTION IS GREATER THAN 50%, THEN USE THAT RELATIVE 
C2A-----REDUCTION AS THE DAMPING VALUE FOR THE NEXT ITERATION.  THE IDEA 
C2A-----HERE IS WHEN WE ARE IN A LINEAR PART OF THE SOLUTION, THE "RRED" 
C2A-----WILL MOVE TOWARDS VALUES NEAR OR > 1, AND LITTLE OR NO DAMPING 
C2A-----WILL BE APPLIED, WHICH IS APPROPRIATE.
      ELSE
        RRED=(RSQ2-RSQ1)/RSQ2/DDAMP
        IF(RRED .GT. 0.5)THEN
          DDAMP=RRED 
          IF(DDAMP .LT. DLOW) DDAMP=DLOW
C2
C2B-----THE RELATIVE REDUCTION IS LESS THAN 50%, INDICATING THE 
C2B-----SOLUTION IS NOT PROGRESSING ADEQUATELY.  BASED ON HOW FAR FROM
C2B-----"RRED" IS FROM 50%, THE DAMPING IS MORE OR LESS AGGRESSIVE.  
C2B-----NOTE: FOR THIS CASE, THE MAXIMUM VALUE OF DAMP WILL BE 0.3 + DLOW
C2B-----SO SETTING DLOW=0.2 IS APPROPIATE SO DAMP WON'T EXCEED 0.5
        ELSE
          DDAMP=0.075/(0.75-RRED) + DLOW 
        ENDIF
C3
C3-----CHECK TO SEE IF SOLUTION IS OSCILLATING BY CHECKING THE 
C3-----RELATIVE DIFFERENCE BETWEEN THE CURRENT RESIDUAL AND THE 
C3-----RESIDUAL FROM TWO ITERATIONS AGO.  IF THEY ARE W/ IN 10%, THEN
C3-----THE SOLUTION IS OSCILLATING.  
C3-----MAKE SURE DDAMP DOESN'T EXCEED DUP
C3-----IF THE SOLUTION IS OSCILLATING, AND WE ARE APPLYING A SIMILAR
C3-----VALUE OF DAMP (W/ IN 1%) THAT WE APPLIED PREVIOUSLY, THEN WE ARE 
C3-----CAUGHT IN SOME TYPE OF ATTRACTOR. SO GENERATE A RANDOM VALUE OF 
C3-----DDAMP AND HOPE TO GET OUT OF THE BASIN OF THIS ATTRACTOR.  
C3-----NOTE THAT THIS IS VERY MUCH A DESPARATION EFFORT!
C
        OSTERM=(RSQ3-RSQ1)/RSQ3
        IF(ABS(OSTERM) .LT. 0.10) OSCIL =.TRUE.
        IF(DDAMP .GT. DUP)DDAMP=DUP
        IF(OSCIL .AND. ABS(DDAMP-DAMP3)/DAMP3 .LT. 0.01)THEN
          NRAN=MOD(NRAN*IMULT+ICONST,IPER)
          QRAND=FLOAT(NRAN)/FLOAT(IPER)
          DDAMP=0.0001D0+QRAND*(1.5D0-0.D0)
        ENDIF
C4
C4-----STORE VALUES OF DAMP AND RESIDUALS FOR NEXT ITERATION
        RSQ3=RSQ2
        DAMP3=DAMP
      ENDIF
      RSQ2=RSQ1
      DAMP=DDAMP
C5
      RETURN
      END
